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Abstract 

Background: Our previous model of the non-isometric muscle fatigue that occurs during repetitive functional 
electrical stimulation included models offeree, motion, and fatigue and accounted for applied load but not 
stimulation pulse duration. Our objectives were to: 1) further develop, 2) validate, and 3) present outcome measures 
for a non-isometric fatigue model that can predict the effect of a range of pulse durations on muscle fatigue. 

Methods: A computer-controlled stimulator sent electrical pulses to electrodes on the thighs of 25 able-bodied 
human subjects. Isometric and non-isometric non-fatiguing and fatiguing knee torques and/or angles were 
measured. Pulse duration (170-600 us) was the independent variable. Measurements were divided into parameter 
identification and model validation subsets. 

Results: The fatigue model was simplified by removing two of three non-isometric parameters. The third remained 
a function of other model parameters. Between 66% and 77% of the variability in the angle measurements was 
explained by the new model. 

Conclusion: Muscle fatigue in response to different stimulation pulse durations can be predicted during non- 
isometric repetitive contractions. 

Keywords: Functional electrical stimulation (FES), Non-isometric, Muscle fatigue, Mathematical model, Pulse 
duration 



Background 

Functional Electrical Stimulation (FES) protocols use 
combinations of stimulation parameters (train duration, 
interpulse interval, pulse duration, and pulse amplitude) 
to produce functional movements in individuals with 
paralysis due to stroke or spinal cord injury (SCI). Un- 
like physiologically induced neuromuscular activation, 
FES synchronously activates motor units according to 
their current thresholds relative to the local extracellular 
current which is dependent on the distance from the 
electrodes [1,2]. Consequently, the recruitment of motor 
units is random [3,4] as compared to the order followed 
by the central nervous system, which recruits the 
smaller, fatigue-resistant motor units first and the larger, 



* Correspondence: sumerian@udel.edu 

'Biomedical Engineering Program, University of California, Davis, CA 95616, 
USA 

Full list of author information is available at the end of the article 



more fatigable motor units last. When the motor units 
are activated synchronously the body cannot derecruit 
motor units as they fatigue and recruit new fresh motor 
units to replace them [5]. This random recruitment 
order together with synchronous activation are thought 
to be two of the major causes for excessive muscle fa- 
tigue during FES. 

A mathematical model capable of predicting angular 
excursion, angular velocity, and joint torque during fa- 
tiguing contractions as a function of the stimulation 
parameters could be used to mathematically test 
combinations of independent and dependent variables to 
identify stimulation strategies that minimize fatigue. In 
addition, a validated model could predict force during 
fatiguing contractions in situations where force cannot 
be measured easily, such as during general non-isometric 
leg extensions. The term non-isometric indicates that the 
joint angle and thus the length of the musculo-tendon 
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unit continually changes as the muscle contracts and 
relaxes. The phrase general non-isometric indicates that 
the leg is free to move solely in response to muscle forces. 
Although many models of non-isometric non-fatiguing 
contractions [6-8] and isometric fatiguing contractions 
[9-14] have been developed, only two models of non- 
isometric fatiguing contractions in humans appear in the 
literature [15,16]. The model by Marion and colleagues 
[16] is the only one that has been experimentally validated 
to predict non-isometric fatigue in response to electrical 
stimulation. 

In our previous study [16] we were interested in 
predicting non-isometric fatigue when the tension per 
activated motor unit was increased through the applica- 
tion of external loads. A similar situation may occur, for 
instance, in the spinal cord injured population when the 
relative resistive torque at the knee as compared to the 
number of activated motor units in the quadriceps 
increases as atrophy progresses. We are now interested 
in determining whether our non-isometric fatigue model 
can predict angular excursion, angular velocity, and joint 
torque due to stimulation of the quadriceps muscles at 
different pulse durations. This interest stems from the 
following reasons: 1) previous studies suggest that 
torque output can be predictably controlled and fatigue 
minimized by simultaneously controlling stimulation 
pulse duration and frequency during repetitive electrical 
stimulation [17-19], 2) others have shown the effect of 
pulse frequency on isometric fatigue and suggest that 
frequency should be minimized [20,21], 3) our isometric 
force-fatigue model accounts for pulse frequency and 
pattern [10,20], but neither the isometric nor the non- 
isometric force-fatigue model account for pulse dur- 
ation, 4) studies suggest that relative isometric fatigue 
(compared to the initial torque) does not change with 
pulse duration [4,21], therefore pulse duration can be 
increased to maintain torque, and 5) the relationship be- 
tween pulse duration and non-isometric fatigue has not 
been reported, therefore it is unknown whether pulse 
duration can be increased to maintain torque and/or ex- 
cursion. Because the overall objective of an ideal FES 
pulse train is to obtain the desired force and motion 
while minimizing fatigue, a fatigue model that takes 
pulse duration into account is required. 

The objectives of this study were to: 1) further develop 
our model of FES non-isometric fatigue to take into 
account pulse duration, while simultaneously minimizing 
the number of parameter identification sessions with 
subjects by minimizing the number of model parameters, 
2) experimentally validate the model at different pulse 
durations, and 3) present outcome measures, such as 
predicted angular excursion, angular velocity, joint torque, 
and power (torque (N-m) x angular velocity (rad/s)) due to 
stimulation, that can be compared over time for different 



independent variables. For consistency with our previous 
study, we chose general non-isometric leg extensions to 
further develop our model of non-isometric muscle fatigue. 
For reference, the term leg is defined as that section of the 
lower limb between the knee and ankle. 

Methods 

Mathematical model 

The force-motion-fatigue model developed by Ding and 
colleagues [7,10,16] was used for this study (see Table 1 
for definitions of symbols). The force-motion model 
[7,16] describes muscle activation, contraction dynamics, 
the force-angle relationship, and the force-angular velocity 
relationship. The input is the time the pulses are delivered, 
and the output is the force (F) at the ankle predicted for 
each time point (see Appendix). 
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Equation 1 models the rate-limiting step that leads to 
the formation of strongly bound crossbridges, and it 
represents the activation dynamics. Equation 2 describes 
the generation of the instantaneous force (F) near the 
ankle due to stimulation. It was derived from a Maxwell 
model of linear viscoelasticity in series with a motor 
[22]. The terms A and G represent the torque-angle [23] 
and torque-angular velocity [8] relationships, respect- 
ively. A torque-pulse duration relationship has not been 
derived yet for this force-motion model of non-isometric 
leg extensions. To meet the overall objective of the 
current study, to predict the effect of a range of pulse 
durations on muscle fatigue, the initial non-fatigue 
torque was measured at the pulse duration of interest 
just prior to the fatigue test. The Michaelis-Menten 
term, CV {K m + Cm), scaled by A and G, drives the devel- 
opment of force. The last term in Equation (2) accounts 
for the force decay over two time constants, and t 2 - 
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Table 1 Definition of symbols and acronyms 


Term 


Unit 


Definition 


Ago 


N ms"' 


Scaling factor reflecting magnitude offeree at 90° 


a 


deg" 2 


Defines parabolic shape of ankle force - knee angle relationship 


0-A 


ms" 2 


Force scaling factor in fatigue model for force-motion model parameter A 


ClKm 


ms^N" 1 


Force scaling factor in fatigue model for force-motion model parameter K m 


Orl 


N" 1 


Force scaling factor in fatigue model for force-motion model parameter Ti 


b 


deg 1 


Defines parabolic shape of ankle force-knee angle relationship 


Pa 


ms" deg" 


Angular velocity x force scaling factor in fatigue model for force-motion model parameter A 




deg"V 


Angular velocity x force scaling factor in fatigue model for force-motion model parameter Km 




ms deg^N" 1 


Angular velocity x force scaling factor in fatigue model for force-motion model parameter ti 


CFT 


_ 


Constant frequency train 


C N 


_ 


Normalized concentration of Ca 2+ -troponin complex 


F 


N 


Instantaneous force near the ankle due to stimulation 


Fload 


N 


Load applied at ankle during general non-isometric leg extensions 


Ffvi 


N 


Represents the resistance to knee extension due to the weight of the leg and all other passive resistance about the knee joint 


TTI 


N s 


Torque Time Integra 




kg m 2 


Net mass moment of inertia of the leg plus the applied load 


Kn 


_ 


Similar to Michaelis-Menten constant. Affinity of actin strong binding site for myosin 


L 


m 


Effective moment arm from knee joint center of rotation to resultant force vector near ankle 


A 


deg 


90° minus the knee flexion angle of the resting non-isometric leg 


n 


_ 


Number of stimuli in train before time t 


Ro 


- 


Characterizes the magnitude of enhancement in C N from the following stimuli 


R, 


_ 


Accounts for differences in activation for each pulse relative to first pulse of train 


SCI 


_ 


Spinal Cord Injury 


ti 


ms 


Time of the I th stimulation 


Tl 


ms 


Time constant offeree decline in the absence of strongly bound cross-bridges 


T 2 


ms 


Time constant offeree decline due to actin-myosin friction in cross-bridges 


Tc 


ms 


Time constant controlling the rise and decay of C N 


Tfat 


ms 


Time constant for force-motion model parameters A, K ml , and T q during fatigue 


e 


deg 


Knee flexion angle, where full extension was 0° 


v, 


N deg" 2 


Scaling factor in the term G 


v 2 


deg" 1 


Constant 


VFT 




Variable frequency train 



Equation 3 models the dynamics of the leg distal to the 
knee. The term F M represents the resistance to knee 
extension due to the weight of the leg and all other 
passive resistance about the knee joint, whereas F load is 
the load applied at the ankle (e.g. 4.54 kg; see Appendix). 
The term A is added to the angle at the knee to ensure 
that angular acceleration is zero at the beginning of 
stimulation. Often the resting knee angle is not exactly 
90°, A is the difference. 

The fatigue model [10,16] monitors changes in the 
three force-motion model parameters that change with 
fatigue, A 90 , K m and Tj. For each time step the input is 



instantaneous force (Equation 2) and angular velocity 
(from angular acceleration in Equation 3) from the 
force-motion model (once all force and fatigue model 
parameters have been identified) for that time step. The 
output is the A 90> K m and r, to be used in the force- 
motion model at the next time step. 



dAgo Ago — A 900 



K m — K m i + K m 2 (5) 
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dK m \ _ K m i — Kmij 
dt 

dK m 2 



Tfat 



dey 

dt ) 



dt 

dr\ 
dt 



-ag m F 



Tl — Tlfl 



Tfat 



(6) 
(7) 

(8) 



The time constant Xf at characterizes the rate of change 
of parameters A 90l K m i, and 7/7 from the pre-fatigue 
values {A 900 , K mli0 , and r, j0 ) to that in a steady state of 
fatigue. All of these terms have been reported previously 
[10,16,24] and the same procedures were used here to 
identify the values. 

Parameter identification 

The force-motion-fatigue model contains a total of 
nineteen parameters. Parameters Rq and t c were held 
constant at 2 (unitless) [10] and 20 ms [20], respect- 
ively (see citations for results showing the derivation of 
these values). Fourteen of the remaining parameters, 
A 90l a, b, Km, Ti, T2, Vi, V2, L/I, and F M: from the force- 
motion model and a A , a Km , a Tl , and Xr at from the 
fatigue model, required identification to both develop 
and validate the model, as well as to generate predictions. 
These parameters were identified from leg extension 
measurements, first from the development then from 
the validation groups of subjects (see Experimental 
Procedures and Figures 1 and 2). The remaining fatigue 
model parameters, Pa, fii<m> and were initially identified 
from measurements and only from the development 
subjects. Model parameters were identified through 
minimization of the sum of squares error between the 
measured and modeled values via a Particle Swarm 
Optimization algorithm [25] followed by a nonlinear 
least-squares algorithm (MatLab") [26]. Optimizations 
were repeated several times to confirm that solutions 
had converged to the "global" minimum. 

Preliminary results showed that for many subjects 
none of the three /3 parameters were needed for accurate 
predictions of the measured angles and angular velocities. 
After careful examination of the subjects that required 
/>, we discovered that only ji Tl was necessary to predict 
fatigue in those subjects. Thus, although ji A and fi Km 
were employed in previous work [16], we postulated 
that j$A and ji Km were not necessary for modeling non- 
isometric fatigue; we explored this hypothesis as described 
in the Results. 

Experimental procedures 
Equipment and participant setup 

Twenty-five healthy subjects, 14 men and 11 women 
(ages 21-48), with no history of lower extremity 



Setup 



Attached electrogonio meter, electrodes, & 
dynamometer or 4.54 kg load to lower limb 



Set pulse amplitude using two 50CFTs, PD 600 us & 170 us 



Stimulated with 14CFTs before each of following tests to potentiate 



Isometric Tests 

Attached to dynamometer; PD = 600 p:s 



Non-Fatigue 

train duty cycle: 1/11 s 



4 knee angles: 15°, 40°, 65°, 90° 

4 trains/angle: 2x50CFT & 2xl2.5VFT 

Model parameters identified @ 90° : A^, K m & %, 

Model parameters identified from all angles: a & b 



Fatigue 

Train duty cycle: 1/2 seconds, except 1st pair 



Knee angle: 90° 
Stim trains: 225 total = 15 sets of lx50CFT, lxl2.5VFT, 

13x33CFT, but 1 st 50CFT & 12.5VFT are non-fatigue 
Model parameters: 14 x (Ago, K„„ & x } ); 0C4, 0* 7W OC^, & Xf u , 



Non- Isometric Tests 

1 train/extension 

+ 



Non-Fatigue 

Train off: 10 seconds 



IsoVelocity 

i 



General - Free Swing 



Velocity: 150 7s 
PD: 600 jis 

Trains: 2x50CFT & 2xl2.5VFT 
Model param: V, & V 2 



Applied load: 4.54 kg 
PD: 170, 200, 250, 400, or 600 jis 
Trains: 2x50CFT & 2xl2.5VFT 
Model param: A w , K m , T,, Z//, & F M 



Fatigue 

5 tests total, 1 test per session 
Train off: 1.2 s, except 1st pair 



Applied load: 4.54 kg 

PD: 170, 200, 250, 400, 600 us; 

Stim trains: 225 total = 15 sets of lx50CFT, lxl2.5VFT, 
13x33CFT, but l" 50CFT & 12.5VFT are non-fatigue 



Figure 1 Flow chart of the testing protocol and the model 
parameters identified from each test. The train duty cycle for the 
general non-isometric contractions varied with pulse duration (PD; 
see results). Parameter r, was identified separately using the force at 
the end of each contraction. Parameters a and b were identified by 
fitting parameter A predicted by Equation (2a) to parameter A from 
Equation (2) for all four knee angles. Fitting Ago, K m , and t, predicted 
by the fatigue model to Ago, K m , and t, from the force model for the 
isometric pre-fatigue and isometric fatiguing contractions identified 
the fatigue model parameters a A , a Km , a lh and j fat (Equations 4-8). 
Initially, parameter (i Tl (Equation 8) was identified by fitting model 
predicted angles and angular velocities to the measurements 
collected at 170, 200, and 600 us during the fatiguing leg 
extensions. An equation for /3 r , (Equation 9) was then derived from 
correlations with the parameters in the model. 
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Force 
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Fatigue 
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1 
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parameters 



Identify fatigue 
parameters 





Force 




Fatigue 




Predicted 




Measured 


B Model 


Model 




Model 




Fatigue 




Fatigue 


Validation f 




1 




1 


i 



3 muscle 
parameters updated 



Validate 
Fatigue Model 



Figure 2 Block diagram demonstrating parameter identification 
(A) and prediction of fatigue (B) using the force-motion (Force) 
and fatigue models. During parameter identification (A) force- 
motion model parameters (muscle parameters) were identified by 
fitting the modeled forces, angles, and angular velocities to the 
measurements collected prior to and during the fatiguing protocol 
in response to the 50CFT and 12.5VFT trains. The fatigue model 
parameters were identified by fitting the parameters A 90 , K m , and r, 
(3 of the 15 muscle parameters) derived from the force-motion 
model to the parameters A go , K m , and r, predicted by the fatigue 
model. During model validation (B), the force and velocity predicted 
by the force-motion model enter the fatigue model at a given time 
step. The fatigue model predicts the parameters Ago, K m , and r, to 
be used by the force-motion model for the next time step. Upon 
completion of all time steps the predicted fatigue is compared to 
the measured fatigue to validate the model. 



orthopedic problems voluntarily participated in this 
study and signed informed consent agreements. This 
study was approved by the University of California 
Human Subjects Review Board. Data from 5 men and 5 
women (ages 19-25) from the previous study on 
predicting fatigue at different loads [16] were also 
analyzed to further validate the model. 

The experimental setup was similar to that described 
previously [16,27] (Figure 1). Subjects were seated in a 
backward-inclined (15° from vertical) chair of an exer- 
cise dynamometer (System 2, Biodex Medical Systems, 
Inc., Shirley, New York). The trunk, hips, and thigh 
were strapped to the chair, thus fixing the hip angle 
and limiting leg movement. The ankle was strapped to 
the lever arm of the dynamometer for the isometric 
and isovelocity tests. The axis of rotation of the knee 
joint was aligned with the axis of rotation of the 
dynamometer. A custom built electrogoniometer with 
two potentiometers, one positioned at the hip and the 
other at the knee axis of rotation, was strapped to the 
lower limb and trunk to measure joint angles. Customized 
software (LabView 8.0, National Instruments Corporation, 
Austin, TX) collected the digitized voltage signals at 300 
Hz from the dynamometer torque transducer and the 
electrogoniometer. Two 7.5 cm x 12.5 cm self-adhesive 
stimulating electrodes (WF35 from www.tensproducts. 
com) were placed on the skin of the right thigh, one at the 



proximal and the other at the distal end of the quadriceps 
muscles. The electrode positions were adjusted until both 
a maximum amplitude and a constant shape of the 
torque-time curve were achieved at 4 different knee 
angles, 90°, 65°, 40° and 20° (where full extension was 0°) 
and at 3 of the 5 pulse durations to be tested (min, mid, 
and max) and until the amount of non-planar movement 
of the leg during general non-isometric leg extensions was 
minimized. In some subjects this resulted in the anode 
being positioned proximal to the cathode. 

Customized software controlled the rate that monophasic 
pulses were delivered by the Grass S48 stimulator 
(Grass Technologies, Astro-Med, Inc. Product Group, 
West Warwick, RI) to the electrodes. A constant-voltage 
transcutaneous system was used to minimize the risk of 
high current densities that can occur with constant- 
current systems if electrode contact with the skin is 
reduced. Others, also studying pulse duration, have used a 
similar system [4,28] . Stimulus efficacy may have changed 
with increasing muscle contraction during delivery of the 
train because the tissues under the skin can move relative 
to the electrodes as the leg extends and because the 
current was not held constant, i.e. maximum stimulation 
of excitable tissue frequently occurs at the beginning of 
the pulse when current is maximum [29,30]. Stimulus 
efficacy also may have changed over time during deliv- 
ery of repetitive fatiguing trains of pulses because of 
sweating, which reduces skin impedance, and because 
of increased blood flow due to increased tissue 
temperature [31]. However, the parameters for the 
force-motion-fatigue model are identified from experi- 
mental measurements from each subject, therefore the 
model can and does account for each subject's muscle 
response to the stimulation system used for the 
measurements. An attached SIU8T stimulus isolation 
unit (Grass Technologies) isolated the electrodes from 
ground, providing greater safety to the subject. 

Testing sessions - general information common to 
All tests 

Each subject participated in 4 to 6 testing sessions. 
Thirteen subjects were used for model development; the 
remaining 12 for model validation. Subjects were asked 
to refrain from strenuous exercise 24 hours prior to each 
testing session. Successive sessions were separated by at 
least 48 hours to allow the muscles to recover and again 
yield the maximum torque measured by the dynamom- 
eter prior to fatigue. Prior to, within the consent to 
participate form, and during the testing sessions, 
participants were asked to relax their legs so that the 
stimulation trains could be applied to relaxed quadri- 
ceps femoris muscles. The consent form states that the 
sessions may have to be repeated if they are unable to 
fully relax their leg. Torque and knee angle were 
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monitored in real-time during the tests. The traces for 
each had consistent shapes appearing at timed intervals 
(during electrical stimulation). Volitional activation 
could be detected easily as alterations to the regularity/ 
uniformity of the traces. Additionally, volitional activa- 
tion during general non-isometric contractions prevents 
or alters the pendulum motion of the leg that occurs 
immediately after the leg drops to the resting position 
after cessation of a stimulation train. If the real-time 
traces on the plots or the pendulum motion of the leg 
looked unusual the test was stopped and the subject was 
gently reminded to relax. 

The stimulation amplitude was set to produce max- 
imum excursion of the freely swinging leg for both the 
minimum and maximum pulse durations while a 4.54 kg 
load was strapped to the ankle. This load was applied 
during all general non-isometric tests. Ten pounds or 
4.54 kg was chosen for two reasons: 1) because it was 
used in previous studies to develop the force-motion 
model used in the current study and to identify its 
parameters [7,16] and 2) because our previous study 
[16] and pilot measurements suggested this load would 
provide measureable declines in force for the desired 
range of pulse durations during the fatiguing contractions. 
The stimulation amplitude was set at the voltage that 
extended the leg to -15° with two 50 Hz trains, one with 
600 us pulses, trains no shorter than 0.2 seconds, and the 
other with 170 us pulses, trains no longer than 0.8 
seconds. This assured a maximum range of motion for 
every subject at all pulse durations, and thus a maximum 
range of fatigue for model development. When using a 
Grass stimulator quadriceps force approaches steady state 
near a pulse duration of 600 us [4,32], therefore 600 us 
was selected as the maximum pulse duration. The mini- 
mum train duration was set at 0.2 seconds so that at least 
two pulses would be delivered at the lowest frequency 
tested. Pulse durations shorter than 170 us were not used 
because the target excursion could not be reached at 
shorter durations by all subjects at the maximum pulse 
amplitude that was limited by the 0.2 second train of 600 
us pulses. Increasing the train duration longer than 0.8 
seconds with 170 us pulses did not increase the excursion 
of the leg. The pulse amplitude depended on the subject, 
ranging from 30 to 83 volts. 

Both constant (CFT) and variable (VFT) frequency 
trains, containing equally spaced singlet pulses or an ini- 
tial doublet (5 ms between pulses within the doublet) 
followed by equally spaced singlet pulses, respectively, 
were applied. Previous studies [10] have shown these 
two types of trains to be effective for identifying the 
model parameters, in particular 50CFT-12.5VFT pairs, 
where 50 and 12.5 refer to the frequency (Hz) of the 
singlet pulses. At the beginning of every test, the quadri- 
ceps were held isometric and stimulated with twelve 



14CFTs (14 Hz pulse frequency), with 0.8 second train 
durations and 5 seconds between trains, to potentiate 
the muscle [23]. Twitch responses initially increase during 
repeated low-frequency stimulation (staircase phenomenon 
or twitch potentiation) and after a tetanic contraction 
(post-tetanic potentiation) [33]. The mechanism of force 
enhancement may be related to phosphorylation of myosin 
light chains and increased Ca 2+ sensitivity [34]. 

Isometric tests 

Non-fatigue isometric Parameters A 90l Km, T lt Tz a, 
and b were identified from the non-fatiguing isometric 
contractions from one testing session. Torque in re- 
sponse to two pairs of testing trains (2 x 50CFT-12 
.5VFT pair) was measured at each of 4 knee angles (15°, 
40°, 65°, 90°). The order of the angles varied from session 
to session and subject to subject. The pulse duration 
was 600 us, train duration was 1 second, and the rest be- 
tween trains was 10 seconds. The muscles rested 4 
minutes between angles, which was sufficient because 
the duty cycle and the number of trains delivered were 
too low to fatigue the muscles [10]. Measured forces were 
compared to modeled forces for initial identification of 
A 90> K m Tj, and t 2 - Parameter A identified from the force- 
motion model (Equation 2) and parameter A predicted by 
the parabolic equation were compared to identify a and b 
(Equation 2a). 

Fatiguing isometric Parameters a a, a-i<m> <*ti< and Tf at 
were identified from the fatiguing isometric contractions 
from one testing session. One fatiguing stimulation 
protocol was applied per subject, at the end of a ran- 
domly selected testing session. The knee angle was 90° 
and the pulse duration was 600 us. Fifteen pairs of 
testing (50CFT-12.5VFT) and 195 fatiguing [33CFT (33 
Hz)] trains, a total of 225 trains were applied as follows: 
1 pair of testing trains followed by 13 fatiguing trains 
and then repeating the 15 trains 15 times. All train 
durations were 1 second. The 50CFT and 12.5VFT in 
the first pair were each followed by a 10 second rest. All 
remaining inter-train rests were 1 second. The 33CFT 
and this duty cycle were chosen because both have been 
proven effective to fatigue the quadriceps within 10 
minutes with minimal discomfort to the participants [10]. 
The 50CFT-12.5VFT pairs, applied after every 13 fatiguing 
trains, generated the forces used for identification of the 
isometric fatigue model parameters [10]. The 15 sets of 
A^ K m and Tj parameters, derived by minimizing the 
error between the forces measured for each pair of testing 
trains and the forces predicted by the force-motion 
model (Equation 2), were compared to the 15 sets of 
A90, K m and Ti parameters predicted by the fatigue 
model (Equations 4-8) to identify a a, a Knv a T j, and Tj at 
(Figure 2). 
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Non-isometric tests 

Non-fatigue Non-isometric Identification of the pa- 
rameters Vj and V 2 in Equation 2b required isovelocity 
measurements from one testing session during which the 
exercise dynamometer extended the leg in passive mode. 
A previous study showed that force-motion model pre- 
dictions were more accurate when parameters Vj and Va 
were identified at 2007second rather than 1257second or 
slower velocities [8]. Therefore, the dynamometer in the 
current study was set to 1507second, its maximum 
velocity in passive mode, and the leg was moved 
from ~110° to 4°. To obtain only the force due to stimula- 
tion, F, it was necessary to collect measurements from leg 
extensions without and with stimulation [8] as the dyna- 
mometer extended the leg from 85° to 20°. This range of 
motion excluded the acceleration and deceleration tails 
and is within the general non-isometric range of motion of 
the leg. Four trains were applied, one per leg extension, two 
50CFT-12.5VFT pairs with pulse durations of 600 us and a 
10 second rest between each train. Measured forces were 
compared to the modeled forces (Equation 2) for identifica- 
tion of Vj, and V 2 . 

Identification of parameters L/I and F M (Equation 3) 
required general non-isometric non-fatiguing measurements 
immediately before every non-isometric fatiguing session. 
The leg was released from the dynamometer, a 4.54 kg load 
was strapped to the ankle, and the leg swung freely. 
Potentiation trains were applied to the free swinging leg im- 
mediately before the general non-isometric measurements. 
Two pairs of testing trains (2 x 50CFT-12.5VFT), each 
followed by a 10 second rest, were applied to the free 
swinging leg, immediately prior to the fatiguing trains in 
the fatigue protocol. The train duration was set to the time 
needed for the leg with attached 4.54 kg load to extend to 
10-15° while the thigh was stimulated with a 50CFT at the 
pulse duration of interest. Measured and modeled angles 
and angular velocities were compared for every non- 
isometric session to identify not only the values for L/I and 
F M (Equation 3), but also to identify the initial, non-fatigue 
force-motion model parameters, A 90 ,a K m i,o, and T2,o> for 
the fatigue model (Equations 4-8), thereby adjusting for 
day-to-day variability. 

Five pulse durations were tested: 170, 200, 250, 400, 
and 600 us, one per testing day. Previous studies [4,32] 
measured the greatest changes in force at pulse durations 
between 100 us and 250 us, at frequencies used in the 
current study. The minimum pulse duration in the current 
study was set to 170 us because shorter pulse durations fre- 
quently did not produce sufficient excursion of the leg at 
the amplitude set for the subject (as described above). The 
next higher pulse duration was set to 200 us because the 
greatest changes in peak force occurred at the lowest pulse 
durations. This small increase in pulse duration produced 
at least a 5% increase in peak force, as was observed in the 



previous studies. The average train durations were 0.64, 
0.51, 0.36, 0.29, and 0.24 seconds for the pulse durations: 
170, 200, 250, 400, and 600 us, respectively. The train 
duration for a given pulse duration was held constant for 
all pulse frequencies. 

Fatiguing Non-isometric Five general non-isometric fa- 
tiguing stimulation protocols were applied per subject, 
one per testing day, immediately following the non-fatigue 
protocol (Figure 1). As with the non-isometric non- 
fatiguing tests, the leg swung freely with a 4.54 kg load 
strapped to the ankle and the same five pulse durations 
were tested: 170, 200, 250, 400, and 600 us. As with the 
isometric fatiguing protocol fifteen pairs of testing 
(50CFT-12.5VFT) and 195 fatiguing [33CFT (33 Hz)] 
trains, a total of 225 trains were applied. The 50CFT 
and 12.5VFT in the first pair were each followed by a 10 
second rest and were used for identification of the initial 
parameters, A 90l o, K m i,o, and t 1i0 , for the fatigue model 
(Equations 4-8) as was stated in the non-fatiguing non- 
isometric section. All remaining inter- train rests were 1.2 
seconds, the minimum time required for the leg to return 
to the resting position (80° to 90°) and to manually stop 
the oscillations with one's hands. The train duration 
remained constant during each fatigue protocol, that is, all 
225 trains for a specific pulse duration test had the same 
train duration. 

Parameters ji A and ji Km were removed from the fatigue 
model and an equation for fi rl (Equation 8) was derived 
during model development from correlations between 
the fitted fi T i and other force-motion-fatigue model 
parameters (Objective 1). Predictions for some subjects 
improved when all three /? parameters were set to 0. 
Therefore, values for fi A , j3 Km , and /> T y were estimated 
separately through optimizations where predictions from 
the fatigue model, containing just one f> per optimization, 
either ji A , fi Km> or /3 T j, were fit to the fatigue measurements 
to determine if one or more /? parameters could be 
eliminated. Preliminary results suggested that ji A and /3 /Cm 
could be removed from the fatigue model. The remaining 
fi T i was initially identified by optimizing the fit between 
the fatigue model values and the angle and angular vel- 
ocity fatigue measurements for the 170, 200, and 600 us 
pulse duration tests. This fitted ji T i was used in the 
correlations to derive an equation for ft Tl . 

Prediction of outcome measures -experimental data from 
both the current and previous study 

Predicted angular excursion, joint torque due to stimula- 
tion, angular velocity, and power (torque (N-m) x angu- 
lar velocity (rad/s)) were compared over time and under 
different pulse duration and load conditions. Two pulse 
durations from the current study and two loads from 
our previous study [16] were used for the comparisons. 
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From the previous study 4.54 kg and 9.08 kg were 
chosen. The 4.54 kg was selected because this was used 
in the current study for all pulse durations and the 9.08 
kg was selected because this was the upper limit. From 
the current study, the pulse durations 600 us and 170 us 
were chosen because these were the lower and upper 
limits tested. A higher pulse amplitude was required in 
the current study than in the previous study to extend 
the leg to -15° at the lowest pulse duration, 170 us. 

Statistical analysis 

To validate the model, the predictive accuracy of the 
model was determined by analysis of the linear regres- 
sion coefficient of determination (r 2 , Objective 2). For 
each subject and each pulse duration in the current 
study (170, 200, 250, 400, and 600 us) or applied load in 
the previous study [16] (0, 1.82, 4.54, 6.36, and 9.08 kg), 
the dependent variable was the predicted, and the 
independent variable was the measured, angular excur- 
sion or angular velocity. Both a fixed slope of unity and 
a y-intercept of zero were used. Ideally, if the predict- 
ive accuracy of the model were 100%, then the linear 
regression r 2 would be unity. Differences in the subject- 
averaged r 2 values between the different pulse durations 
or applied loads, both for angular velocity and excursion 
were determined using repeated measures ANOVAs 
followed by Tukey post hoc tests. A two-factor test was 
used for the subjects tested in the current study where the 
independent variables were pulse duration (170, 200, 250, 
400, and 600 us, non-isometric and isometric) and type of 
subject (development and validation). A one-factor test 
was used for the subjects tested in the previous study 
where the independent variable was load (0, 1.82, 4.54, 
6.36, and 9.08 kg). In all cases the dependent variable was 
the r 2 -value. 

To present outcome predictions (Objective 3), differences 
in predicted angular excursion, torque time integral (TTI), 
joint torque at maximum power, angular velocity at 
maximum power, and maximum power due to stimulation 
of the quadriceps were determined using two-factor 
repeated measures ANOVAs followed by Tukey post 
hoc tests. The independent variables for the two-factor 
ANOVAs were pulse duration (170, 200, 250, 400, and 600 
us; measured in the current study) or load (0, 1.82, 4.54, 
6.36, and 9.08 kg; measured in the previous study [16]) and 
contraction number (the first 33CFT and the average of the 
last seven 33CFTs). The 33 Hz train was chosen because it 
was used to fatigue the muscle and was the middle 
frequency train, between the 50 Hz and 12.5 Hz trains. The 
last seven trains were averaged because the torque-time 
and angle-time curves typically varied more at the end of 
the fatigue protocol than at the beginning. Additionally, at 
the beginning of the fatigue protocol there was a 10 second 
rest just prior to the first 33CFT, whereas only 1.2 seconds 



separated the remaining trains in the fatigue protocol. The 
shorter rest time resulted in somewhat increased variability 
in the starting position and velocity before each contrac- 
tion. Because the fatigue curve was at steady state when the 
last set of 33CFTs was applied, the average of the last half 
of that set adequately represented the last train. The 
dependent variables were predicted angular excursion, TTI, 
joint torque at maximum power, angular velocity at max- 
imum power, and maximum power, all due to stimulation. 
The predicted joint torque was computed by multiplying 
the force predicted by the force-motion-fatigue model by 
the moment arm (L) from the knee joint center of rotation 
to the center of the load applied just proximal to the ankle. 
In all cases p < 0.05 was considered significant. 

Results 

Modifications to the fatigue model (objective 1) 

Complete data sets were collected on 25 subjects. Pre- 
liminary regression analyses of predictions of fatigue 
using the force-motion-fatigue model from the previous 
study [16], which used equations for f} A , fi Knv and (> Tl 
from the fatigue model (Equations 4-8), showed that al- 
though this model accounted for most of the variance in 
most subjects, predictions for some subjects improved 
when all three fi parameters were set to 0 (not shown). 
Preliminary results suggested that inclusion of parameter 
fi Tl alone, without fi A or f> Km , in the fatigue model could 
account for fatigue in all the subjects. Angular velocity 
multiplied by /S rJ (Equation 8) reduced the impact of fa- 
tigue model parameter a Tl on the force relaxation time 
constant Ti. Parameter a rl accounts for the increase in 
Ti that occurs during isometric fatigue, but in some 
subjects, parameter t x changed less during non-isometric 
fatigue than during isometric fatigue (Figure 3 shows an 
extreme case). Applying this new fatigue model to 
measurements from our previous study [16] confirmed 
that Pa and ji Km were not needed in the fatigue model to 
predict non-isometric fatigue. 

The parameter fi Tl could be expressed as a function of 
parameters in the non-isometric force and isometric fa- 
tigue models. This is shown in Equation 9: 

P T1 = 8.5XHT 11 

A 1.6 v p0.5 v T 3.5 v -,1.3 v T 0.9 
^90,0 * r M * n.O.iso x "ti x L fat /n , 

x ' (9) 

where A 90 ,o> F M , and Vj are non-fatigue force-motion 
model parameter values from the day of the non- 
isometric fatigue session of interest, t^ojso is the non- 
fatigue force-motion model parameter value from the 
isometric fatigue session, and a Tl and tf at are fatigue 
model parameter values from the isometric fatigue ses- 
sion. The equation for /? rJ (Equation 9) in the current 
study is different from the equation in the previous 
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Figure 3 Effects of fatigue model parameters a T , and ft, on isometric and non-isometric contractions in one subject where ft 7 was 
higher than average. (A) Addition of ft,(d0/dt) to a Tl in Equation 8 brings force relaxation time constant, ii, closer to pre-fatigue value (top), 
resulting in non-isometric predictions with a faster rate of fatigue (bottom). Thin solid black lines are measured values, just prior to (80°) and at 
end of extension. (B) For isometric contractions, removal of a T , in Equation 8 keeps t q constant at pre-fatigue value, resulting in isometric 
predictions with a faster rate of fatigue (FTI =force time integral). (C) A single isometric contraction shows that when a Tl is included in isometric 
Equation 8, progressively slower twitch relaxation times increase the force of contraction. (D) Single non-isometric contractions show that 
addition of ft,(d0/dt) to Equation 8 partially negates the effect of a Tj . In A. and B. pairs of 50CFT and 1 2.5VFT testing trains were followed by 
13x33CFT fatiguing trains. Contractions in (C) and (D) occurred at 0.7 min (dashed line) and 2.3 min (solid line). Non-isometric: 4.54 kg load, 250 
us pulse duration. Isometric: 600 us pulse duration. Initial force-motion model parameters: Ago = 2.10 N/ms, K m = 3.52e-01, r, = 36.1 ms, i 2 = 52.1 
ms, t c = 20 ms, R 0 = 2, a = -4.49e-004 deg~ 2 , b = 3.44e-02 deg" 1 , V, = 3.71e-01 N/deg 2 , V 2 = 2.29e-02 deg -1 , LA = 9.85 kg'W 1 , F M = 247.5 N. 
Fatigue model parameters: r fat = 99.4 s, a A = -4.03e-07 ms" 2 , a Km = -1.36e-08 ms~V, a T , = 2.93e-05 N'\ ft, = 8.54e-04 ms deg^KF 1 . 



study (Equation 8a) (11) because in the previous study 
three /? parameters, pA, p Km , and ji Tl , were used in the 
fatigue model. All three were identified simultaneously 
when fitting the fatigue model (Equations 4-8) predictions 
to the fatigue measurements to obtain the fitted ft values 
used in the correlations to derive the equations for ji. In 
the current study only one ji parameter, /3 rJ (Equation 8), 
was used and identified when fitting the fatigue model 
predictions to the measurements, therefore the fitted 
fi Tl in the current study was different from that in the 
previous study. Because li zl could be estimated from 
equation 9, non-isometric fatigue measurements were 
not needed to predict non-isometric fatigue. 

Predictions of fatigue validated the model (objective 2) 

Both measured and predicted angular velocity and 
excursion showed the greatest fatigue at the highest 
load, shortest pulse duration, and longest train duration 
(Figure 4, Objective 2). Train duration was a confounding 
factor, but was consistent across both studies. Predicted 



velocity- and excursion-time curves were within one stand- 
ard deviation of measured curves, with the exception of the 
first 1.5 minutes of the 0 kg load tests (Figure 4B, D). 

Comparison of predictions to measurements through 
linear regression analyses (Figure 5) indicated that the new 
non-isometric force-motion-fatigue model accounted for 
between 66% and 77% of the variability in non-isometric 
fatigue for different clinically relevant pulse durations 
(170, 200, 250, 400, or 600 us) with 4.54 kg applied to the 
ankle (Figure 6A). Predictions of measurements from our 
previous study [16] indicated that the new model also 
explained between 67% and 81% of the variability in 
non-isometric fatigue for different applied loads (0, 1.82, 
4.54, 6.36, or 9.08 kg) when stimulating with 600 us 
pulses (Figure 6B). Recall that the model development 
measurements were collected only in the current study and 
only at 170, 200, and 600 us. All other measurements were 
used only for model validation. The predictions for the iso- 
metric measurements exceeded those for the non-isometric 
measurements (0.0001<p< 0.02, Figure 6B), accounting for 
>85% of the variability in isometric fatigue. 



Marion et al. Journal of NeuroEngineering and Rehabilitation 2013, 10:13 
http://www.jneuroengrehab.eom/content/10/1/13 



Page 10 of 16 





01234567 012345678 
Time (min) Time (mln) 



Figure 4 Measured (Ms) and predicted (Pr) pulse duration (A,C) and load (B,D) dependent reduction in relative angular velocity (A and 
B±SD) and excursion (C±SD and D) during fatiguing contractions. Angular excursion is defined as the difference between the initial and final 
knee angle of a leg extension. The 33CFT contractions shown are normalized to the first contraction. Measurements for (B) and (D) were 
collected in our previous study [1 6]. Predictions are within one standard deviation of measurements, with the exception of the first 1 .5 minutes of 
the 0 kg load. In (A) and (C) applied load is 4.54 kg and average train durations in (A), from shortest to longest pulse duration, are 0.64, 0.51, 0.36, 
0.29, and 0.24 seconds. In (B) and (D) pulse duration is 600 us and average train durations in (D), from highest load to lowest load, are: 0.89, 0.54, 
0.51, 0.32, 0.19 seconds. Maintaining a constant excursion necessitated changing the train duration. Only 2 loads and pulse durations are shown 
in (B) and (C), respectively, so that the standard deviations could be shown. 



Outcome measures that can be predicted and compared 
(objective 3) 

Torque at the knee due to stimulation of the quadriceps 
cannot be measured directly during general non- 
isometric leg extensions because the leg is not attached 
to any device that might resist its natural motion. 
However, this torque can be predicted by our force- 
motion-fatigue model. In this way, angular excursion, 
joint torque, angular velocity, and power due to stimula- 
tion can be compared over time and under different 
conditions (Figure 7, Objective 3). The predicted dependent 
variables showed significant fatigue (contraction number as 
the independent variable) at both loads and both pulse 
durations. With applied load or pulse duration as the 



independent variable, differences between the two 
applied loads or two pulse durations were not always 
significant. The predicted initial maximum power was 
not significantly different between the two loads or between 
the two pulse durations. The predicted angular velocity at 
maximum power was significantly less at the highest load 
and lowest pulse duration, while the predicted initial joint 
torque at maximum power was significantly greater at the 
highest load and lowest pulse duration. The initial angular 
excursion at 170 us pulse duration was significantly less 
than at 600 us (Figure 7A). Keep in mind that the train 
duration was set such that the 50CFT, not necessarily the 
33CFT, produced the maximum excursion at each pulse 
duration or load. 
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Figure 5 Predicted (Pr) vs. measured (Ms) angular excursion 
from one subject for the 600 us and 170 us pulse durations. 

Applied load was 4.54 kg. Pairs of 50CFT and 12.5VFT testing trains 
were followed by 13 x 33CFT fatiguing trains (15 sets of 15 
contractions for a total of 225 contractions). Row 2 is the linear 
regression analysis. Both a fixed slope of unity and y-intercept of 0 
were used, because this is the ideal relationship between 
measurement and prediction. Initial force-motion model parameters 
for the 600 us and 1 70 us pulse durations, respectively, were: A go = 
1.11 and 1.69 N/ms, K m = 2.75e-01 and 3.60e-01, r, = 57.9 and 31.5 
ms, t 2 = 59.8 ms, J c = 20 ms, R 0 = 2, a = -3.27e-04 deg~ 2 , b = 4.21 e- 
02 deg" 1 , V, = 1 .56 N/deg 2 , V 2 = 4.98e-02 deg" 1 , Ul = 22.2 and 5.86 
kg" m" , Fm = 86.3 and 254.7 N. Fatigue model parameters were: 
T fat = 95.4 s, a A = -4.02e-07 ms" 2 , a Km = -7.34e-08 ms"V, 
a T , = 4.17e-05 N" 1 , ft, = 4.60e-05 and 1.53e-04 ms deg~V\ 



Discussion 

The key findings in the current study were that 

(a) pulse duration was not explicitly needed in the 
fatigue model; its effects on fatigue were captured 
by its effects on force, 

(b) two of three /3 parameters could be eliminated from 
our previous fatigue model without loss of 
predictive value with current and previous data sets, 

(c) the remaining ji parameter is expressed completely 
as a function of values already measured, so, 
effectively, no additional parameters were added to 
the fatigue model, 

(d) the new force-motion-fatigue model accounted for 
66-77% and 67-81% of the variability in the non- 
isometric measurements from the current and 
previous study, respectively, and 

(e) the model can be used to compare the power, angular 
velocity, angular excursion, and joint torque due to 
stimulation produced during fatiguing non-isometric 
contractions under different testing conditions. 




600 400 250 200 170 0 1.8 4.5 6.4 9.1 
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Figure 6 Average linear regression coefficients of 
determination (r 2 ; ± 95% confidence limit) for predicted versus 
measured angular excursion and velocity of all contractions 
(50CFT, 12.5VFT and 33CFT). The non-isometric force-motion- 
fatigue model accounted for 66-81% of the variability in fatigue 
during general non-isometric leg extensions. Isometric force-time 
integral at 90° (iso) is shown for comparison. In A. n of gray bars = 
13 model development subjects (note that only 170, 200, and 600 
us were used for model development) and n of black bars = 12 
model validation subjects. Applied load was 4.54 kg. In B. n = 10 
model validation subjects (pulse duration = 600 us). % - compared 
to 0 kg (0.001 <p<0.05); - compared to all 5 pulse durations 
(0.0001 <p< 0.02). Because the greatest potentiation occurred during 
the 0 kg load tests and because the force-motion-fatigue model 
does not include a term for potentiation, predicted excursions and 
velocities were lower than the measured values for 0 kg. 



The fatigue model was simplified by eliminating the 
parameters yS A and f> Km from the fatigue model and gen- 
erating a new equation for (3 rl (Equation 9) as a function 
of existing force-motion-fatigue model parameters. Be- 
cause f$ Tl was multiplied by negative angular velocity, 
the yS T i term reduced the effect of a T j, bringing Tj closer 
to its pre-fatigue value (Figure 3). In some subjects the 
difference between the pre-fatigue and fatigue twitch re- 
laxation times was minimal during non-isometric 
contractions. For some subjects, the twitch relaxation 
time during non-isometric fatiguing contractions was 
less than during isometric fatiguing contractions. 

Non-isometric fatigue measurements were not needed 
to predict non-isometric fatigue. In total, all but 5 
parameters {A 90 , K m T lt L/I and F M ), from both the force 
and fatigue models were identified from measurements 
collected during one testing session. The remaining 5 
parameters were identified from pre-fatigue general non- 
isometric leg extension measurements from each non- 
isometric fatigue testing session. 

The predictive ability of our new non-isometric force- 
motion-fatigue model (0.66 <= r <= 0.77 for pulse duration 
and 0.67 <= r 2 <= 0.81 for applied load) tended to be higher 
than that of our previous non-isometric force-motion-fa- 
tigue model (0.56 < r <= 0.76 for applied load) [16], though 
lower than that of our isometric (r 2 >0.85) force-fatigue 
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Figure 7 Predicted angular excursion, torque time integral 
(TTI), joint torque at maximum power, angular velocity at 
maximum power, and maximum power, all due to stimulation 
(mean ± SD; n = 25 (A) and 10 (B) subjects; 33CFT). Note that 
initial maximum power was not significantly different between the 
two loads or between the two pulse durations, but velocity was 
lowest at the highest load and shortest pulse duration. Isometric TTI 
(iso; 90°) is shown for comparison. Gray bars = first 33CFT; black bars 
= average of last half of last set of 33CFTs. Average train durations 
for 600 us, 170 us, 4.5 kg, 9.1 kg, and isometric (600 us pulse 
duration) were: 0.24, 0.64, 0.51, 0.89, and 1.00 seconds, respectively. 
e - compared to first contraction (p<0.0001); ip q - compared to 4.5 
kg and 9.1 kg (0.0001 < p<0.02), or to 600 us and 170 us (p<0.0001); 
upo - compared to 4.5 kg (0.0001 < p< 0.03) or to 600 us (0.0001 < 
p< 0.02). 



model (Figure 6). The predictions in the current study for 
the measurements collected in the previous non-isometric 
modeling study [16] tended to be more accurate than those 
in the previous study because 1) the baseline angle or 
velocity for every stimulation train (contraction) delivered 
during the fatigue protocol on a given day in the current 
study was the initial value before the first train in the 
fatigue protocol, whereas in the previous study the baseline 
was an average of the initial angles or velocities before each 
of the 225 trains, and these sometimes deviated from the 
baseline of the resting leg and 2) the time between the last 
potentiation train and the first train in the fatigue protocol 
was reduced in the current study compared to the previous 
study, which reduced the magnitude of force enhancement 
that often occurred within the first few trains in the 
fatigue protocol. Insufficient potentiation explains why 



the measured fatigue tended to be less than the 
predicted fatigue for the 0 kg load (Figure 4) because 
the force-motion-fatigue model had no provision for 
potentiation. 

A number of factors may explain why the isometric 
force-fatigue model accounted for more of the variability 
in the isometric measurements (Figure 6; 86-92%) than 
the non-isometric force-motion-fatigue model could ac- 
count for in the non-isometric measurements. These in- 
clude the following: 

(1) In the isometric case, all model parameters were 
identified from force measurements at one knee 
angle, 90°. In the non-isometric case, the force- 
length relationship model parameters were 
identified from force measurements at 4 different 
angles and the isovelocity and free model 
parameters were identified from angle and angular 
velocity measurements at the angles between ~85° 
(resting) and -12° (nearly full extension). 

(2) In the isometric case, the electrode position relative 
to the nerves and muscles beneath was nearly 
constant from the beginning to the end of a 
fatiguing protocol. In the non-isometric case, both 
the skin and muscles moved as the leg extended, 
and maximum extension depended on pulse 
frequency and extent of fatigue, and therefore the 
amount of movement may have varied from train to 
train. 

(3) In the isometric case, the leg remained in the 
sagittal plane. In the non-isometric case the leg may 
have moved out of the sagittal plane as it fatigued. 

(4) In the isometric case, the potentiation protocol 
given just prior to the fatigue protocol minimized 
the force enhancement that often occurred during 
the first few trains in the fatigue protocol. In the 
non-isometric case, the potentiation protocol was 
not as effective at reducing the force enhancement 
that occurred with the shortest train durations and 
highest velocities, perhaps due to the continual and 
rapid change in myofiber or myofibril conformation. 

(5) In the isometric case, the initial force immediately 
before every contraction in the fatigue protocol was 
the same. In the non-isometric case, the initial angle 
and angular velocity was not always identical 
because we manually stopped and released the leg 
after each fatiguing extension, allowing for some 
human error. 

To reach the desired excursion, as pulse duration 
decreased, train duration increased; as applied load 
increased, train duration increased. Train duration was 
therefore a confounding factor in our results, interacting 
with pulse duration and applied load (Figure 4). Taken 
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together, the results of both studies suggest that the 
higher the duty cycle of the train, the greater the fatigue 
(constant rest time between trains: 1.2 and 1.3 seconds). 
This has been observed by others [35]. It may seem that 
holding the train duration the same across all pulse 
durations (or loads) would have led to a clearer inter- 
pretation of the measurements, but then maximum ex- 
cursion would not have been constant across trials. Both 
the 12.5 Hz VFT and 50 Hz CFT trains were required to 
identify the pre-fatigue force-motion model parameters. 
Considering that the leg was free to move, the maximum 
train duration was limited by the highest frequency train 
at the longest pulse duration. Holding the train durations 
constant would have resulted in significantly different 
angular excursions among pulse durations, thus creating a 
different confounding factor. Additionally, our objective 
was to validate fatigue predictions from excursion and 
velocity measurements. Using the same pre-fatigue 
excursion (at 50 Hz) for every pulse duration provided 
the largest range of excursion between the pre-fatigue 
and final fatigue measurements. 

Comparing initial and final outcome measures in re- 
sponse to different independent variables, such as applied 
load or pulse duration, could help a therapist determine 
which stimulation parameters are most desirable for the 
patient and task. If higher joint torque is required (e.g. to 
strengthen the muscles), then a pulse duration of 170 us is 
preferable to 600 us (see Figure 7). On the other hand, if 
maintaining the highest level of power over the greatest 
length of time is the goal, then a pulse duration of 600 us 
is preferable to 170 us. There was no significant difference 
in the initial maximum power between the two pulse 
durations however, the final maximum power for the 170 
us was less than that for 600 us. 

The isometric force model has been shown to perform 
equally well for both able-bodied and SCI subjects, 
requiring only minor modifications to the parameter 
identification procedures for the SCI subjects [28,36]. 
The new non-isometric force-motion-fatigue model, also 
validated to account for different loads per activated 
muscle (which could occur if atrophy progresses), may 
be equally robust, where similar minor modifications to 
the parameter identification procedures would pertain to 
this model. The maximum force generating ability of the 
muscles could be estimated from peak twitch force 
measurements as described by Ding, et al. (2005) [36]. 
The stimulation amplitude could be set as described in 
the current study, but would not exceed a level consist- 
ent with 50% of the force generating ability of the 
muscles. The isometric experimental protocol and iden- 
tification of the isometric force model parameters could 
be similar to that described by Ding, et al. (2005) [36]. 
The non-isometric experimental protocol could be simi- 
lar to that described in the current study, but the pulse 



frequencies would be as described by Ding, et al. (2005) 
[36]. The model parameters are subject specific, identified 
by fitting the model to the experimental measurements 
obtained from one testing session; therefore the current 
procedure for identifying these model parameters may 
require only minor modifications for the non-isometric 
force-motion-fatigue model to predict fatigue in SCI 
subjects. Spastic measurements would be excluded. 

Our model has the potential to help physical therapists 
design stimulation protocols for patients in rehabilitation 
programs and to help researchers improve the task per- 
formance of FES systems [19,32,37-39]. The isometric 
force-fatigue model was extensively validated to account 
for the effect of different pulse frequencies and patterns 
on fatigue [10,20]. The non-isometric force-motion-fa- 
tigue model has been validated to account for different 
applied loads and pulse durations, and these have 
resulted in a validation of different train duty cycles. 
From these model validations we learned that frequency, 
pulse pattern, pulse duration, and applied load are not 
explicitly needed in the fatigue model. Their effects on 
fatigue can be captured by their effects on force. There- 
fore, the non-isometric force-motion-fatigue model 
should be able to predict unique combinations of 
stimulation parameters for different subjects, such that 
each subject can achieve a desired outcome, such as 
maintaining a functional level of power for a useful 
period of time (e.g. Figure 7). The non-isometric force- 
motion model [39] and the isometric fatigue model 
[40] have been used in a similar manner in other studies. 
The force-motion-fatigue model, with all model parameters 
identified for the task, could either mathematically test 
combinations of stimulation parameters until the desired 
outcome is obtained, or could be fit to an experimental 
force or trajectory (using an optimization algorithm) to 
generate optimal stimulation patterns that yield the force or 
trajectory for the desired length of time (see Maladen, 
et al. [39]). 

Because this non-isometric force-motion-fatigue 
model would be capable of generating subject-specific 
and task-specific stimulation patterns that can maintain 
a desired force and motion for a desired length of time 
into the future, it has the potential for use as a feed for- 
ward model in FES systems [41]. If a system either does 
not use a feed forward model or requires more immedi- 
ate real time output, then this model could be used to 
test the performance of the system prior to patient use. 
The model could generate a series of task-specific opti- 
mal stimulation patterns, and these patterns could be 
compared to the real time FES system selections to 
optimize the system. 

Because able-bodied subjects were tested in our study, 
there is a small chance that volitional contractions oc- 
curred during stimulation. However, it is unlikely that 
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volitional contractions, if present, substantively affected 
our results. The force-motion-fatigue model has been 
shown to successfully predict fatigue in response to differ- 
ent frequencies and pulse patterns for numerous subjects 
over many years [10,27,36,40]. This indicates that the 
signal-to-noise ratio has been high, where here the 
stimulated contractions correspond to the signal and the 
volitional contractions correspond to the noise. In all 
cases, several testing sessions were performed on each 
subject, and each session was separated by 48 hours. 

Conclusion 

Pulse duration was not explicitly needed in the fatigue 
model; its effects on fatigue were captured by its effects 
on force. The non-isometric force-motion-fatigue model 
from our previous study [16] was simplified to predict 
non-isometric fatigue both at different applied loads and 
at different pulse durations. Parameters p A and fi Km in 
the previous version of the fatigue model were eliminated 
and a new equation for the parameter f> Tl was derived. The 
($ Tl was solely a function of existing model parameters; 
therefore measurements of non-isometric fatigue are not 
needed to predict non-isometric fatigue. From 66% to 
77% of the variability in the non-isometric measurements 
for different pulse durations was explained by the new 
force-motion-fatigue model. This new non-isometric 
force-motion-fatigue model can be used to predict angular 
excursion, angular velocity, joint torque, or power due to 
stimulation at different time intervals during repetitive 
contractions. This could assist with rehabilitation 
exercises and with the design and testing of new FES 
control systems. 

Appendix 

Derivation of the equation of motion 

As described by Perumal, et al [8] the instantaneous mo- 
ment about the knee center of rotation was derived from 
the free body diagram of the leg shown in Figure 8. 

The equation of motion derived from the free body 
diagram for isovelocity extensions is: 



FbwL 



-mg kos(8)-H+ T st 



(10) 



where F ext = F Bio is the force component of the torque 
measured by the Biodex dynamometer and 

Tstim — FL 



Thus, 



(10a) 



(10b) 



where F is the force just proximal to the malleoli exerted 
by the quadriceps through the knee joint in response to 
stimulation. It is defined as the instantaneous force near 



Lstim 




Figure 8 Schematic representation of the leg, modeled as a rigid 
body. L is the distance from the center of the knee joint to either the 
center of the calf pad of the Biodex dynamometer knee attachment or 
the center of the load applied just proximal to the ankle when the leg 
is not attached to the dynamometer (just proximal to the malleoli but 
distal to the prominent calf musculature), / is the distance from the 
center of the knee joint to the center of mass of the tibia, J stim is the 
torque at the knee due to stimulation, F ext is the either the force 
measured by the Biodex dynamometer {F Bio ) if the leg is attached to 
the dynamometer (the resistance that the force dynamometer exerts 
against the ankle to maintain a constant angular velocity) or the 
component of the applied ankle weight (F^od) that resists the 
contractile force of the quadriceps if the leg is swinging free, mg is the 
weight of the leg below the knee and the foot, and H is the resistance 
moment to knee extension due to the visco-elasticity of the structures 
at the knee. 



the ankle due to stimulation in Table 1. Previous passive 
force measurements on healthy subjects showed that 
(see Perumal, et al [8]): 



^ = 7?cos(<9) 

where R is an intermediate variable. 
Letting 



(10c) 



mgl 



R 



and substituting equation lOd into 10b yields 

F = F Bio + F M cos(0) 



(lOd) 



(lOe) 



where F M is obtained by fitting the function F M cos(6)) to 
force data collected during passive leg extensions where 
the quadriceps are relaxed and the dynamometer alone 
extends the leg. 
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The equation of motion for the general non-isometric 
leg extensions is: 

= F load Lcos(6) + mglcos(9)+H- FL (11) 

~ = j{(F[oad+F M )cos(e)-F} (11a) 

Angular acceleration is no longer zero. F ext = the com- 
ponent of Fi oad (applied ankle weight) that resists the 
contractile force of the quadriceps. The parameter L/I is 
a lumped parameter encompassing more than length 
and moment of inertia. Previous estimates of L/I using 
anthropometric data revealed differences from the values 
estimated through optimization [7,8]. The differences 
may be the result of: 1) identifying L/I and F M simultan- 
eously during optimization and 2) assuming that acceler- 
ation and/or applied weight have no effect on the force- 
motion model (Equation 2), keeping in mind that F in 
the equation of motion is predicted from the force- 
motion model (Equation 2). Therefore, L/I represents a 
more generalized parameter. 
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